---
title: "A Comment on 'Instrumentally Inclusive: The Political Psychology of Homonationalism' (@turnbull2024instrumentally)."
author: 
  - name: "Daniel de Kadt"
    affiliation: "London School of Economics and Political Science"
    email: "d.n.de-kadt@lse.ac.uk"
citeproc: true
bibliography: references.bib
format: html


---

## Setup

Note to replicators: As much as possible, I have attempted throughout to avoid adjusting the original authors' code. Most changes are documented as comments. Please note I have used `groundhog` -- please run the setup chunk once to install the relevant packages with the correct versions. 

```{r setup, echo = FALSE, include = FALSE, warnings = FALSE}
# load libraries
library(groundhog)

pkgs <- c("tidyverse", "jtools", "ggpubr", "ggrepel", "patchwork", "gt", "modelsummary", "interactions", "margins", "skimr", "survey", "estimatr")
groundhog::groundhog.library(pkgs, "2024-06-01")

# create output directory if it doesn't exist
if(!dir.exists("output")) dir.create("output")

# modelsummary options
options(modelsummary_model_labels = "Computer Modern")

# set seed exactly as per replication materials
set.seed(1)
```

```{r data ingestion, echo = FALSE, include = FALSE, warnings = FALSE}
# load both datasets
uk <- read_csv("study1_data.csv")
load("study2_data.Rda")

# set color palette as per replication materials
colors<- c("#205C8A", "#d11141")

# cleaning per replication materials
uk <- uk%>% 
  mutate(treat= as.factor(treatment),
         treatnum= as.numeric(treatment),
         gender= as.factor(gender),
         degree= as.factor(degree),
         nonwhite= as.factor(nonwhite),
         queer= as.factor(queer),
         relig= as.factor(relig),
         religion= as.factor(religion),
         race= as.factor(race),
         fourarm= as.factor(fourarm),
         immbelow= as.factor(immbelow),
         imm3= as.factor(imm3),
         region= as.factor(region),
         voterecall= as.factor(voterecall),
         brexit= as.factor(brexit),
         ideology= as.factor(ideology),
         agecat= as.factor(agecat))

# cleaning per replication materials
spain <- spain%>% 
  mutate(treat= as.factor(treat),
         treatnum= as.numeric(treat),
         gender= as.factor(gender),
         supportcat= as.factor(support),
         agecat= as.factor(agecat),
         child= as.factor(child),
         immdum= as.factor(immdum),
         imm5= as.factor(imm5),
         imm3= as.factor(imm3),
         foreignborn= as.factor(foreignborn),
         CCAA= as.factor(CCAA),
         queer= as.factor(queer))
```

```{r uk analysis and subsetting, echo = FALSE, include = FALSE, warnings = FALSE}
# subset uk data per replication materials
treat <- subset(uk, treatnum==1)
control <- subset(uk, treatnum==0)
proimm <- subset(uk, immbelow==0)
noproimm <- subset(uk, immbelow==1)

# analyse uk data per replication materials
modelsub1<- lm(support ~ treat, data=proimm)
proimm$predictedb <- predict(modelsub1, proimm)

modelsub2<- lm(support ~ treat, data=noproimm)
noproimm$predictedb <- predict(modelsub2, noproimm)

# subset again per replication materials
treatsub1 <- subset(proimm, treatnum==1)
controlsub1 <- subset(proimm, treatnum==0)
treatsub2 <- subset(noproimm, treatnum==1)
controlsub2 <- subset(noproimm, treatnum==0)
```

```{r spain analysis and subsetting, echo = FALSE, include = FALSE, warnings = FALSE}
# subset spain data per replication materials
treatES <- subset(spain, treat==1)
controlES <- subset(spain, treat==0)
proimmES <- subset(spain, immdum==1)
noproimmES <- subset(spain, immdum==0)

# analyse spain data, but use lm() not glm() so that Spain and UK analyses are exactly the same:
modelsub1ES <- lm(support ~ treat, weight=nationalweight, data=proimmES)
proimmES$predictedb <- predict(modelsub1ES, proimmES)

modelsub2ES <- lm(support ~ treat, weight=nationalweight, data=noproimmES)
noproimmES$predictedb <- predict(modelsub2ES, noproimmES)

# subset again per replication materials
treatsub1ES <- subset(proimmES, treat==1)
controlsub1ES <- subset(proimmES, treat==0)
treatsub2ES <- subset(noproimmES, treat==1)
controlsub2ES <- subset(noproimmES, treat==0)
```

## Figure 1: Weight Distribution in Study 2 (Spain)

```{r spain weights histogram, echo = FALSE, include = FALSE, warnings = FALSE}
# create weights plot
weights_plot <- ggplot(spain, aes(x=nationalweight)) +
  geom_histogram(bins=30, fill="black", color="black") +
  theme_minimal() +
  labs(x = "Weight",
       y = "Frequency") 

```

```{r spain weights histogram plot, echo = FALSE, warnings = FALSE}
#| label: fig-weights
#| fig-cap: "Distribution of Weights in Study 2  (Spain)"

ggsave(weights_plot, filename = "output/weights_plot.jpg", width=7, height=5)

weights_plot

```

## Table 1: Summary of Variance Estimation Throughout Paper and Supplementary Materials

```{r se-table, echo = FALSE}
#| label: tbl-se
#| tbl-cap: "Summary of variance estimation throughout paper and supplementary materials."

df <- data.frame(
  Study = c("Study 1:", rep("",6),
            "Study 2:", rep("",7)
            ),
  Item = c(
  "", "Figure 3", "Figure 4", "Figure A4", "Table A7", "Table A8",
  "",
  "", "Figure 6", "Figure 7", "Figure 8", "Figure A5", "Table A9", "Table A10", "Table A11"
  ),
  RobustSEs = c("", "No", "Yes", "No", "No", "Yes",
                "",
                "", "No", "No", "No", "No", "No", "No", "Yes"),
  Notes = c(
    "", #study 1:
    "90% CIs shown",
    "Text overlay uses classical SE, 90% CIs shown",
    "90% CIs shown",
    "", "", 
    "",
    "", #study 2:
    "90% CIs shown",
    "Text overlay uses classical SE, 90% CIs shown",
    "90% CIs shown",
    "90% CIs shown",
    "", "", ""  
  ),
  stringsAsFactors = FALSE
)

# Create the gt table

tbl_se <- df %>%
  gt() %>%
  cols_align(
    align = "left",
    columns = c(Study, Item, Notes)
  ) %>%
  cols_label(
    Study = "",
    Item = "Figure/Table",
    RobustSEs = "Robust SEs?",
    Notes = "Notes"
  ) %>%
  tab_options(
    table.font.size = "small"  
  ) 

# Save as .tex file:
latex_tbl_se <- as_latex(tbl_se)
writeLines(as.character(latex_tbl_se), "output/tbl_se.tex")

# Show the table
tbl_se

```

## Table 2: Reproduction of Table A9 and Summary of Researcher Choices

```{r spain replication no-weights, echo = FALSE, include = FALSE}
# computational reproduction of Table A9:
models_rep <- list(
  'Base' = lm(support ~ treat + imm_1, data=spain, weight = nationalweight),
  'Interaction model' = lm(support ~ treat*imm_1, data=spain, weight = nationalweight),
  'Pro-immigration only' = lm(support ~ treat, data=proimmES, weight = nationalweight),
  'Anti-immigration only' = lm(support ~ treat, data=noproimmES, weight = nationalweight)
)

# reproduction by analysis, with varying researcher choices:
models_base <- list(
  'Replication' = lm(support ~ treat + imm_1, data=spain, weight = nationalweight),
  'No Weights' = lm(support ~ treat + imm_1, data=spain),
  'Weighted Robust SE' = estimatr::lm_robust(support ~ treat + imm_1, data=spain, weight = nationalweight, se_type = "HC3"),
  'Unweighted Robust SE' = estimatr::lm_robust(support ~ treat + imm_1, data=spain, se_type = "HC3"),
  'Weighted Survey-Robust SE' = survey::svyglm(support ~ treat + imm_1, design=svydesign(ids=~1, weights=~nationalweight, data=spain))
)

models_int <- list(
  'Replication'  = lm(support ~ treat*imm_1, data=spain, weight = nationalweight),
  'No Weights' = lm(support ~ treat*imm_1, data=spain),
  'Weighted Robust SE' = estimatr::lm_robust(support ~ treat*imm_1, data=spain, weight = nationalweight, se_type = "HC3"),
  'Unweighted Robust SE' = estimatr::lm_robust(support ~ treat*imm_1, data=spain, se_type = "HC3"),
  'Weighted Survey-Robust SE' = survey::svyglm(support ~ treat*imm_1, design=svydesign(ids=~1, weights=~nationalweight, data=spain))
)

models_proimmes <- list(
  'Replication'  = lm(support ~ treat, data=proimmES, weight = nationalweight),
  'No Weights' =  lm(support ~ treat, data=proimmES),
  'Weighted Robust SE' = estimatr::lm_robust(support ~ treat, data=proimmES, weight = nationalweight, se_type = "HC3"),
  'Unweighted Robust SE' = estimatr::lm_robust(support ~ treat, data=proimmES, se_type = "HC3"),
  'Weighted Survey-Robust SE' = survey::svyglm(support ~ treat, design=svydesign(ids=~1, weights=~nationalweight, data=proimmES))
)

models_noproimmes <- list(
  'Replication'  = lm(support ~ treat, data=noproimmES, weight = nationalweight),
  'No Weights' =  lm(support ~ treat, data=noproimmES),
  'Weighted Robust SE' = estimatr::lm_robust(support ~ treat, data=noproimmES, weight = nationalweight, se_type = "HC3"),
  'Unweighted Robust SE' = estimatr::lm_robust(support ~ treat, data=noproimmES, se_type = "HC3"),
  'Weighted Survey-Robust SE' = survey::svyglm(support ~ treat, design=svydesign(ids=~1, weights=~nationalweight, data=noproimmES))
)

# ancillary mechanism test (A11):
mech <- list(
  'EU norms' = estimatr::lm_robust(pride_valoresUE ~ treat*imm_1, data=spain, se_type = "HC3"),
  'Western liberal values' = estimatr::lm_robust(pride_libertadOCC ~ treat*imm_1, data=spain, se_type = "HC3"),
  'Green politics' = estimatr::lm_robust(pride_verde ~ treat*imm_1, data=spain, se_type = "HC3"),
  'Domestic violence protections' = estimatr::lm_robust(pride_viomach ~ treat*imm_1, data=spain, se_type = "HC3"),
  'Spanish flag' = estimatr::lm_robust(pride_bandera ~ treat*imm_1, data=spain, se_type = "HC3"),
  'Spanish military efforts' = estimatr::lm_robust(pride_mili ~ treat*imm_1, data=spain, se_type = "HC3")
)

```


```{r table-summary, echo=FALSE}
#| label: tbl-summary
#| tbl-cap: "Summary of Researcher Choices, Point Estimates, and Statistical Conclusions. "

get_stars <- function(est, se) {
  t <- abs(est / se)
  stars <- if (t > 2.576) "***"
           else if (t > 1.96) "**"
           else if (t > 1.645) "*"
           else ""
  c(paste0(round(est, 3), stars), paste0("(", round(se, 3), ")"))
}

get_stars_interaction <- function(est1, se1, est2, se2) {
  stars1 <- if (abs(est1/se1) > 2.576) "***"
            else if (abs(est1/se1) > 1.96) "**"
            else if (abs(est1/se1) > 1.645) "*"
            else ""
  stars2 <- if (abs(est2/se2) > 2.576) "***"
            else if (abs(est2/se2) > 1.96) "**"
            else if (abs(est2/se2) > 1.645) "*"
            else ""
  est_str <- paste0(round(est1, 3), stars1, " $\\mid$ ", round(est2, 3), stars2)
  se_str <- paste0("(", round(se1, 3), ") $\\mid$ (", round(se2, 3), ")")
  c(est_str, se_str)
}

table_data <- data.frame(
  Choices = c(
    "Weighted + Classical SE", 
   "",
    "Unweighted + Classical SE", 
   "",
    "Weighted + Robust SE",  
   "",
    "Unweighted + Robust SE", 
   "",
    "Weighted + Survey-R SE",
   ""
  ),
  
  Baseline = c(
    get_stars(summary(models_base$Replication)$coefficients[2,1], summary(models_base$Replication)$coefficients[2,2]),
    get_stars(summary(models_base$`No Weights`)$coefficients[2,1], summary(models_base$`No Weights`)$coefficients[2,2]),
    get_stars(summary(models_base$`Weighted Robust SE`)$coefficients[2,1], summary(models_base$`Weighted Robust SE`)$coefficients[2,2]),
    get_stars(summary(models_base$`Unweighted Robust SE`)$coefficients[2,1], summary(models_base$`Unweighted Robust SE`)$coefficients[2,2]),
    get_stars(summary(models_base$`Weighted Survey-Robust SE`)$coefficients[2,1], summary(models_base$`Weighted Survey-Robust SE`)$coefficients[2,2])
  ),
  
  Interaction = c(
    get_stars_interaction(
      summary(models_int$Replication)$coefficients[2,1], summary(models_int$Replication)$coefficients[2,2],
      summary(models_int$Replication)$coefficients[4,1], summary(models_int$Replication)$coefficients[4,2]),
    get_stars_interaction(
      summary(models_int$`No Weights`)$coefficients[2,1], summary(models_int$`No Weights`)$coefficients[2,2],
      summary(models_int$`No Weights`)$coefficients[4,1], summary(models_int$`No Weights`)$coefficients[4,2]),
    get_stars_interaction(
      summary(models_int$`Weighted Robust SE`)$coefficients[2,1], summary(models_int$`Weighted Robust SE`)$coefficients[2,2],
      summary(models_int$`Weighted Robust SE`)$coefficients[4,1], summary(models_int$`Weighted Robust SE`)$coefficients[4,2]),
    get_stars_interaction(
      summary(models_int$`Unweighted Robust SE`)$coefficients[2,1], summary(models_int$`Unweighted Robust SE`)$coefficients[2,2],
      summary(models_int$`Unweighted Robust SE`)$coefficients[4,1], summary(models_int$`Unweighted Robust SE`)$coefficients[4,2]),
    get_stars_interaction(
      summary(models_int$`Weighted Survey-Robust SE`)$coefficients[2,1], summary(models_int$`Weighted Survey-Robust SE`)$coefficients[2,2],
      summary(models_int$`Weighted Survey-Robust SE`)$coefficients[4,1], summary(models_int$`Weighted Survey-Robust SE`)$coefficients[4,2])
  ),
  
  ProImmigration = c(
    get_stars(summary(models_proimmes$Replication)$coefficients[2,1], summary(models_proimmes$Replication)$coefficients[2,2]),
    get_stars(summary(models_proimmes$`No Weights`)$coefficients[2,1], summary(models_proimmes$`No Weights`)$coefficients[2,2]),
    get_stars(summary(models_proimmes$`Weighted Robust SE`)$coefficients[2,1], summary(models_proimmes$`Weighted Robust SE`)$coefficients[2,2]),
    get_stars(summary(models_proimmes$`Unweighted Robust SE`)$coefficients[2,1], summary(models_proimmes$`Unweighted Robust SE`)$coefficients[2,2]),
    get_stars(summary(models_proimmes$`Weighted Survey-Robust SE`)$coefficients[2,1], summary(models_proimmes$`Weighted Survey-Robust SE`)$coefficients[2,2])
  ),
  
  AntiImmigration = c(
    get_stars(summary(models_noproimmes$Replication)$coefficients[2,1], summary(models_noproimmes$Replication)$coefficients[2,2]),
    get_stars(summary(models_noproimmes$`No Weights`)$coefficients[2,1], summary(models_noproimmes$`No Weights`)$coefficients[2,2]),
    get_stars(summary(models_noproimmes$`Weighted Robust SE`)$coefficients[2,1], summary(models_noproimmes$`Weighted Robust SE`)$coefficients[2,2]),
    get_stars(summary(models_noproimmes$`Unweighted Robust SE`)$coefficients[2,1], summary(models_noproimmes$`Unweighted Robust SE`)$coefficients[2,2]),
    get_stars(summary(models_noproimmes$`Weighted Survey-Robust SE`)$coefficients[2,1], summary(models_noproimmes$`Weighted Survey-Robust SE`)$coefficients[2,2])
  )
)

# If there are any factors, convert them to character first to avoid issues
table_data_final <- table_data
table_data_final[] <- lapply(table_data_final, as.character)

# Create the table
tbl_summary <- table_data_final |>
  gt() |>
  cols_align(
    align = "center",
    columns = c(Baseline, Interaction, ProImmigration, AntiImmigration)
  )  %>%
  cols_label(
    Choices = "Researcher Choices",
    Baseline = "Base Model",
    Interaction = "Interaction",
    ProImmigration = "Pro-immig.",
    AntiImmigration = "Anti-immig.",
  ) %>%
 tab_style(
  style = cell_borders(
    sides = "top",
    color = "black",
    weight = px(1.5)
  ),
  locations = cells_body(
    rows = nrow(table_data_final) 
    )
  ) %>%
  tab_options(
    table.font.size = "small"  
  ) %>%
  fmt_markdown(
    columns = everything()
  )

# Save as .tex file:
latex_tbl_summary <- as_latex(tbl_summary)
writeLines(as.character(latex_tbl_summary), "output/tbl_summary.tex")

# Show the table
tbl_summary

```

## Tables 3 and 4: 
### 3: Study 2 (Spain): Heterogeneous Effects By Weight Bin for Anti-Immigrant Sub-Sample
### 4: Study 2 (Spain): Heterogeneous Effects By Weight Bin for Pro-Immigrant Sub-Sample

```{r weight bin prep, echo = FALSE, include = FALSE, warnings = FALSE}
# create bins of weights:
spain <- spain %>%
  mutate(weightbins = case_when(
    nationalweight >= 0 & nationalweight < 0.01 ~ "low",
    nationalweight >= 0.1 & nationalweight < 3 ~ "medium",
    nationalweight >= 3 ~ "high",
    TRUE ~ "other"
  ))

# re-subset
proimmES <- subset(spain, immdum==1)
noproimmES <- subset(spain, immdum==0)

# build tables - anti-immigration subsample
models_hetfx_anti <- list(
  'Low Weights' = estimatr::lm_robust(support ~ treat, data=noproimmES[noproimmES$weightbins=="low",], se_type = "HC3"),
  'Mid Weights' = estimatr::lm_robust(support ~ treat, data=noproimmES[noproimmES$weightbins=="medium",], se_type = "HC3"),
  'High Weights' = estimatr::lm_robust(support ~ treat, data=noproimmES[noproimmES$weightbins=="high",], se_type = "HC3")
)

# build tables - pro-immigration subsample
models_hetfx_pro <- list(
  'Low Weights' = estimatr::lm_robust(support ~ treat, data=proimmES[proimmES$weightbins=="low",], se_type = "HC3"),
  'Mid Weights' = estimatr::lm_robust(support ~ treat, data=proimmES[proimmES$weightbins=="medium",], se_type = "HC3"),
  'High Weights' = estimatr::lm_robust(support ~ treat, data=proimmES[proimmES$weightbins=="high",], se_type = "HC3")
)


# build tables - interaction hetfx
models_hetfx_int <- list(
  'Low Weights' = estimatr::lm_robust(support ~ treat*imm_1, data=spain[spain$weightbins=="low",], se_type = "HC3"),
  'Mid Weights' = estimatr::lm_robust(support ~ treat*imm_1, data=spain[spain$weightbins=="medium",], se_type = "HC3"),
  'High Weights' = estimatr::lm_robust(support ~ treat*imm_1, data=spain[spain$weightbins=="high",], se_type = "HC3")
)
```


```{r weight bin hetfx table anti, echo = FALSE, warnings = FALSE}
#| label: tbl-weight-hetfx-anti
#| tbl-cap: "Study 2 (Spain): Heterogeneous Effects By Weight Bin for Anti-Immigrant Sub-Sample"

tbl_weight_hetfx_anti <- suppressWarnings(
  modelsummary(
    models_hetfx_anti, 
    output = "gt", 
    coef_map = c('treat1' = 'Treatment', '(Intercept)' = 'Intercept'),
    gof_omit = "BIC|AIC|R2 Adj.|F|RMSE|Log.Lik.",  
    stars = c('*'=.1, "**"=.05, "***"=.01), 
    add_rows = tribble(
      ~term,  ~Base,  ~Base, ~Base,
      "Weighted", "No", "No", "No",  
      "HC3 Robust SEs", "Yes", "Yes", "Yes",
      "Survey Robust SEs", "No", "No", "No"
    )
  )
  ) %>%
    tab_options(
    table.font.size = "small"  
  ) %>%
  fmt_markdown(
    columns = everything()
    )

# Save as .tex code:
latex_tbl_weight_hetfx_anti <- as_latex(tbl_weight_hetfx_anti)
latex_tbl_weight_hetfx_anti <- gsub(
  "(?s)\\\\begin\\{minipage\\}\\{.*?\\}.*?\\\\end\\{minipage\\}",
  "",
  latex_tbl_weight_hetfx_anti,
  perl = TRUE
)

writeLines(as.character(latex_tbl_weight_hetfx_anti), "output/tbl_weight_hetfx_anti.tex")

# Show the table
tbl_weight_hetfx_anti

```

```{r weight bin hetfx table pro, echo = FALSE, warnings = FALSE}
#| label: tbl-weight-hetfx-pro
#| tbl-cap: "Study 2 (Spain): Heterogeneous Effects By Weight Bin for Pro-Immigrant Sub-Sample"
tbl_weight_hetfx_pro <- suppressWarnings(
  modelsummary(models_hetfx_pro, output = "gt", 
             coef_map = c('treat1' = 'Treatment',  
                          '(Intercept)' = 'Intercept'),
             gof_omit = "BIC|AIC|R2 Adj.|F|RMSE|Log.Lik.",  
             stars = c('*'=.1, "**"=.05, "***"=.01), 
             add_rows = tribble(~term,  ~Base,  ~Base, ~Base,
               "Weighted", "No", "No", "No",  
               "HC3 Robust SEs", "Yes", "Yes", "Yes",
               "Survey Robust SEs", "No", "No", "No", 
             ))
  ) %>%
    tab_options(
    table.font.size = "small"  
  ) %>%
  fmt_markdown(
    columns = everything()
    )

# Save as .tex code:
latex_tbl_weight_hetfx_pro <- as_latex(tbl_weight_hetfx_pro)
latex_tbl_weight_hetfx_pro <- gsub(
  "(?s)\\\\begin\\{minipage\\}\\{.*?\\}.*?\\\\end\\{minipage\\}",
  "",
  latex_tbl_weight_hetfx_pro,
  perl = TRUE
)

writeLines(as.character(latex_tbl_weight_hetfx_pro), "output/tbl_weight_hetfx_pro.tex")

# Show the table
tbl_weight_hetfx_pro
```

## Table 5: Study 2 (Spain): Covariate Balance by Weight Bin

```{r weight balance analysis, echo = FALSE, include = FALSE, warnings = FALSE}
# create a balance table where we show the mean of imm_1, age, edu, gender, and other covariates by weightbin:
bal_table <- spain %>%
  mutate(gender = as.numeric(as.character(gender)) - 1) %>%
  group_by(weightbins) %>%
  summarise_at(
    c("age", "gender", "edu", "child", "foreignborn", "queer", "imm_1"),
    ~round(mean(as.numeric(as.character(.)), na.rm = T),2)
  ) %>%
# transpose the table so the weightbins are the columns, and the covariates the rows:
  pivot_longer(cols = c(age, gender, edu, child, foreignborn, queer, imm_1)) %>%
  pivot_wider(names_from = weightbins, values_from = value) %>%
  relocate(low, .before = high) %>%
  rename(`Low Weight Bin` = low, `Medium Weight Bin` = medium, `High Weight Bin` = high, `Variable` = name)
```

```{r weight balance table, echo = FALSE, warnings = FALSE}
#| label: tbl-weight-balance
#| tbl-cap: "Study 2 (Spain): Covariates by Weight Bin"

tbl_weight_balance <- gt(bal_table) %>%
    tab_options(
    table.font.size = "small"  
  ) %>%
  fmt_markdown(
    columns = everything()
    )

# Save as .tex code:
latex_tbl_weight_balance <- as_latex(tbl_weight_balance)
writeLines(as.character(latex_tbl_weight_balance), "output/tbl_weight_balance.tex")

# Show the table
tbl_weight_balance

```

## Figure 2: Study 1 (UK): Corrected Reproduction of Figure 3

```{r uk interaction analysis, echo = FALSE, include = FALSE, warnings = FALSE}
# correct the code to use a linear regression lm() and not logistic regression. 
model1 <- lm(support ~ treat*imm_1, data=uk)
# use lm_robust for robust SEs that can be used by margins(). Set se_type = "HC3" to be consistent with summ(.,robust=TRUE). Point estimates are of course numerically identical bar rounding.
model1_robust <- estimatr::lm_robust(support ~ treat*imm_1, data=uk, se_type = "HC3")

# create plots per the replication materials, adding in the CI and making SE robust. Must use lm class object here not lm_robust object.
pred<- interact_plot(model1, pred = imm_1, modx = treat, interval = TRUE, robust = TRUE,
                     colors = colors) +
  labs(title="",
       y="Predicted support for\nLGBT+ education in schools",
       x="")+
  theme_minimal()+
  theme(legend.position = "none",
        axis.text.x =element_blank())+
  annotate(
    geom="text", x = 2.5, y = .65, size = 4, color = "#d11141", fontface=2,
    label = "Slope for\ntreated group")+
  annotate(
    geom = "curve", x =1.6, y = .6, xend = 1.2, yend = .44, 
    curvature = .4, arrow = arrow(length = unit(2, "mm")), colour="#d11141")+
  annotate(
    geom="text", x = 6, y = .4, size = 4, color = "#205C8A", fontface=2,
    label = "Slope for\ncontrol group")+
  annotate(
    geom = "curve", x =5, y = .4, xend = 2.52, yend =.38, 
    curvature = -.4, arrow = arrow(length = unit(2, "mm")), colour="#205C8A")

gg_df <-
  # update to robust object
  model1_robust %>%
  margins(at = list(imm_1 = seq(0, 10, by = 1))) %>%
  summary %>%
  as.data.frame() %>%
  filter(factor == "treat1")

ame<- ggplot(gg_df, aes(imm_1, AME)) +
  geom_point(colour="#d11141") +
  geom_line(colour="#d11141") +
  coord_cartesian(xlim = c(0, 10), ylim = c(-.25, .25)) +
  # change to 95% ci:
  geom_errorbar(aes(ymax = (AME-SE*1.96), ymin = (AME+SE*1.96)), width = 0, colour="#d11141") +
  geom_hline(yintercept = 0, linetype = "dashed", colour="#205C8A") +
  xlab("Pre-treatment attitudes towards immigration")+
  ylab("Conditional ATE") +
  theme_minimal()
```

```{r print uk interaction plots, echo = FALSE, warnings = FALSE}
#| label: fig-3-correction
#| fig-cap: "Study 1 (UK): Corrected Reproduction of Figure 3"

fig_3_correction <- pred/ame+ 
  plot_annotation(title = 'Conditional average treatment effect: Study 1 (UK)',
                  theme = theme(plot.title = element_text(size = 14, face="bold")))

ggsave(fig_3_correction, filename = "output/fig_3_correction.jpg", width=8, height=5)

fig_3_correction

```

## Figure 3: Study 1 (UK): Study 2 (Spain): Corrected Reproduction of Figure 6 (Weights Used)

```{r spain interaction analysis, echo = FALSE, include = FALSE, warnings = FALSE}
# correct the code to use a linear regression lm() and not logistic regression. 
modelES <- lm(support ~ treat*imm_1, data=spain, weight=nationalweight)
# use lm_robust for robust SEs that can be used by margins(). Set se_type = "HC3" to be consistent with summ(.,robust=TRUE). Point estimates are of course numerically identical bar rounding.
modelES_robust <- estimatr::lm_robust(support ~ treat*imm_1, data=spain, weight=nationalweight, se_type = "HC3")

# create plots per the replication materials, adding in the CI and making SE robust:
predES<- interact_plot(modelES, pred = imm_1, modx = treat, interval = TRUE, robust = TRUE,
                        colors = colors)+
  labs(title="",
       y="Predicted support for\nLGBT+ education in schools",
       x="")+
  theme_minimal()+
  theme(legend.position = "none",
        axis.text.x =element_blank())+
  annotate(
    geom="text", x = 2.5, y = .65, size = 4, color = "#d11141", fontface=2,
    label = "Slope for\ntreated group")+
  annotate(
    geom = "curve", x =1.6, y = .6, xend = 1.2, yend = .44, 
    curvature = .4, arrow = arrow(length = unit(2, "mm")), colour="#d11141")+
  annotate(
    geom="text", x = 6, y = .4, size = 4, color = "#205C8A", fontface=2,
    label = "Slope for\ncontrol group")+
  annotate(
    geom = "curve", x =5, y = .4, xend = 2.6, yend =.31, 
    curvature = -.4, arrow = arrow(length = unit(2, "mm")), colour="#205C8A")

gg_df <-
  # update to robust object
  modelES_robust %>%
  margins(at = list(imm_1 = seq(0, 10, by = 1))) %>%
  summary %>%
  as.data.frame() %>%
  filter(factor == "treat1")

ameES<- ggplot(gg_df, aes(imm_1, AME)) +
  geom_point(colour="#d11141") +
  geom_line(colour="#d11141") +
  coord_cartesian(xlim = c(0, 10)) +
  # change to 95% ci:
  geom_errorbar(aes(ymax = (AME-SE*1.96), ymin = (AME+SE*1.96)), width = 0, colour="#d11141") +
  geom_hline(yintercept = 0, linetype = "dashed", colour="#205C8A") +
  xlab("Pre-treatment attitudes towards immigration")+
  ylab("Conditional ATE") +
  theme_minimal()
```

```{r print spain interactions plots, echo = FALSE, warnings = FALSE}
#| label: fig-6-correction
#| fig-cap: "Study 2 (Spain): Corrected Reproduction of Figure 6 (Weights Used)"

fig_6_correction <- predES/ameES+ 
  plot_annotation(title = 'Conditional average treatment effect: Study 2 (Spain)',
                  theme = theme(plot.title = element_text(size = 14, face="bold")))

ggsave(fig_6_correction, filename = "output/fig_6_correction.jpg", width=8, height=5)

fig_6_correction

```

## Figure 4: Study 1 (UK): Corrected Reproduction of Figure 4

```{r uk dot-plot processing, echo = FALSE, include = FALSE, warnings = FALSE}
proimmplot<- effect_plot(model = modelsub1, pred = treat, robust=TRUE, 
                         cat.geom="point", cat.interval.geom="linerange",
                         # correct the 95% ci:
                         colors="black", cat.pred.point.size=3, int.width = .95)+
  labs(title = paste0("Pro-immigrant voters (N = ", nrow(proimm),")"))+
  ylab("Support for LGBT+ education in schools (0-1)")+
  xlab("")+
  scale_y_continuous(limits=c(-0.1,1.1), breaks=c(0, .2, .4, .6, .8, 1)) +
  scale_x_discrete(labels=c("0" = "Control", "1" = "Treatment"))+
  # correct the y variable:
  geom_jitter(data=treatsub1, aes(x=treat, y=support),
              height=.05, width=.2, alpha=.35, shape=20,
              pch=21, size=2, color="#d11141")+
  # correct the y variable:
  geom_jitter(data=controlsub1, aes(x=treat, y=support),
              height=.05, width=.2, alpha=.35, shape=20,
              pch=21, size=2, color="#205C8A")+  
  geom_bracket(xmin = c("0"), xmax = c("1"),
               y.position = c(.45), label = c("ATE=-.06*"),
               tip.length =-0.05,
               color="black")+
  theme_minimal()+
  theme(axis.title.y = element_text(face="bold"),
        axis.text.x = element_text(face="bold"))

noproimmplot<- effect_plot(model = modelsub2, pred = treat, robust=TRUE, 
                           cat.geom="point", cat.interval.geom="linerange",
                           # correct the 95% ci:
                           colors="black", cat.pred.point.size=3, int.width = .95)+
  labs(title = paste0("Anti-immigrant voters (N = ", nrow(noproimm),")"))+
  ylab("Support for LGBT+ education in schools (0-1)")+
  xlab("")+
  scale_y_continuous(limits=c(-0.1,1.1), breaks=c(0, .2, .4, .6, .8, 1)) +
  scale_x_discrete(labels=c("0" = "Control", "1" = "Treatment"))+
  # correct the y variable:
  geom_jitter(data=treatsub2, aes(x=treat, y=support),
              height=.05, width=.2, alpha=.35, shape=20,
              pch=21, size=2, color="#d11141")+
  # correct the y variable:
  geom_jitter(data=controlsub2, aes(x=treat, y=support),
              height=.05, width=.2, alpha=.35, shape=20,
              pch=21, size=2, color="#205C8A")+
  geom_bracket(xmin = c("0"), xmax = c("1"),
               y.position = c(.79), label = c("ATE=.10**"),
               tip.length =0.05,
               color="black")+
  theme_minimal()+
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.text.x = element_text(face="bold"))

```

```{r print uk dot-plots, echo = FALSE, warnings = FALSE}
#| label: fig-4-correction
#| fig-cap: "Study 1 (UK): Corrected Reproduction of Figure 4"

fig_4_correction <- proimmplot+noproimmplot+ 
  plot_annotation(title = 'Effect of out-group treatment among:',
                  caption = "Lines are 95% confidence intervals.Treatment group outcome statistically distinct at p<.1(*), p<0.05(**), & p<0.01(***)",
                  theme = theme(plot.title = element_text(size = 14, face="bold")))

ggsave(fig_4_correction, filename = "output/fig_4_correction.jpg", width=8, height=5)

fig_4_correction

```

## Figure 5: Study 2 (Spain): Corrected Reproduction of Figure 7  (Weights Used)

```{r spain dot-plot processing, echo = FALSE, include = FALSE, warnings = FALSE}
proimmplotES<- effect_plot(model = modelsub1ES, pred = treat,
                            cat.geom="point", cat.interval.geom="linerange",
                           # correct 95% ci and make them robust:
                            colors="black", cat.pred.point.size=2, int.width = .95, robust = TRUE)+
  labs(title = paste0("Pro-immigrant voters (N = ", nrow(proimmES),")"))+
  ylab("Support for LGBT+ education in schools (0-1)")+
  xlab("")+
  scale_y_continuous(limits=c(-0.1,1.1), breaks=c(0, .2, .4, .6, .8, 1)) +
  scale_x_discrete(labels=c("0" = "Control", "1" = "Treatment"))+
  # correct the y variable:
  geom_jitter(data=treatsub1ES, aes(x=treat, y=support),
              height=.05, width=.2, alpha=.35, shape=20,
              pch=21, size=2,  color="#d11141")+
  # correct the y variable:
  geom_jitter(data=controlsub1ES, aes(x=treat, y=support),
              height=.05, width=.2, alpha=.35, shape=20,
              pch=21, size=2, color="#205C8A")+  
  geom_bracket(xmin = c("0"), xmax = c("1"),
               # correct significance from *** to ** due to robust SEs:
               y.position = c(.45), label = c("ATE=.11**"),
               tip.length =-0.05,
               color="black")+
  theme_minimal() +
  theme(axis.title.y = element_text(face="bold"),
        axis.text.x = element_text(face="bold"))

noproimmplotES<- effect_plot(model = modelsub2ES, pred = treat, 
                              cat.geom="point", cat.interval.geom="linerange",
                             # correct 95% ci and make them robust:
                              colors="black", cat.pred.point.size=2, int.width = .95, robust = TRUE)+
  labs(title = paste0("Anti-immigrant voters (N = ", nrow(noproimmES),")"))+
  ylab("Support for LGBT+ education in schools (0-1)")+
  xlab("")+
  scale_y_continuous(limits=c(-0.1,1.1), breaks=c(0, .2, .4, .6, .8, 1)) +
  scale_x_discrete(labels=c("0" = "Control", "1" = "Treatment"))+
  # correct the y variable:
  geom_jitter(data=treatsub2ES, aes(x=treat, y=support),
              height=.05, width=.2, alpha=.35, shape=20,
              pch=21, size=2, color="#d11141")+
  # correct the y variable:
  geom_jitter(data=controlsub2ES, aes(x=treat, y=support),
              height=.05, width=.2, alpha=.35, shape=20,
              pch=21, size=2, color="#205C8A")+
  geom_bracket(xmin = c("0"), xmax = c("1"),
               # correct significance from *** to [] due to robust SEs:
               y.position = c(.79), label = c("ATE=.10"),
               tip.length = 0.05,
               color="black")+
  theme_minimal()+
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.text.x = element_text(face="bold"))

```

```{r print spain dot-plots, echo = FALSE, warnings = FALSE}
#| label: fig-7-correction
#| fig-cap: "Study 2 (Spain): Corrected Reproduction of Figure 7  (Weights Used)"

fig_7_correction <- proimmplotES+noproimmplotES+ 
  plot_annotation(title = 'Effect of out-group treatment among:',
                  caption="Lines are 95% confidence intervals. Treatment group outcome statistically distinct at p<.1(*), p<0.05(**), & p<0.01(***)",
                  theme = theme(plot.title = element_text(size = 14, face="bold")))

ggsave(fig_7_correction, filename = "output/fig_7_correction.jpg", width=8, height=5)

fig_7_correction

```
